
FUNCTION ldss3_define_slitlets,exposure_file,norder=norder,width=width

;+----------------------------------------------------------------------------------------
;
;  LDSS3_DEFINE_SLITLETS
;
;  Trace slit edges in 2D spectra to create a 2D array of slit positions.
;  The slit position at a pixel is the x-intercept of the curve of
;  constant wavelength that would pass through the center of that
;  pixel.
;
;  Pretty much finished. Still some tweaking to do. The initial identification of 
;  the peaks (slit edges) could be a problem
;  
;
;
;  INPUT : 
;   - exposure file : flat field image. Typically dome flats
;                     
;  OUTPUT : 
;
;   the main output is slitgrid.fits which gives the position of the 
;   y-intercept (x=0) for any pixel x,y
;
;
;  KEYWORDS : 
;
;  - norder : order of the polynomial fit to the slit edge along the x-direction. 
;             (norder is 3 by default). 
;
;  - width : one-sided width where to search for slit edge from the initial y-location 
;            of the peaks. 
;
;
;  HISTORY : 
; 
;  Written by JRG 10/27/11
;
;
;+-------------------------------------------------------------------------------------------
if (n_params() eq 0) then begin
   print,'CALLING SEQUENCE : result=ldss3_define_slitlets(exposure_file,norder=norder,width)'
   return,-1
endif

; read parameter file : 
@params.list

dat=readfits(exposure_file+'.fits',fhead)
sizef=size(dat)
szx=sizef(1)
szy=sizef(2)
dat=dat*1.0
; display image
;window,0,xsize=szx,ysize=szy
; display the map : 
;tv,dat,TRUE=interleaving[0]
;nb=1

imy=shift(dat,0,-1)-dat 

; create a stack :
nbs=9
imystack=imy(xmax-nbs:xmax,*)
sumstack=abs(total(imystack,1))
idy=indgen(n_elements(sumstack))
derivy=abs(deriv(idy,sumstack))
nbpeaks=0
peaksid=peaks(sumstack,100,npk)
nbpeaks=n_elements(peaksid)
print,'nbpeaks =',n_elements(peaksid)
print,'peaksid =',peaksid

window,1,xsize=900,ysize=400
plot,idy,sumstack,xrange=[min(peaksid),max(peaksid)],yrange=[0,0.9*max(sumstack)]
oplot,idy(peaksid),sumstack(peaksid),psym=5
; oplot,idy,derivy,linestyle=2

nslit=0.5*nbpeaks
print,'Nb of slits =',nslit
slit=replicate({slitstrct},nslit)

ids_eve=indgen(23)*2
lledge=peaksid(ids_eve)
ids_odd=2*indgen(23)+1
uuedge=peaksid(ids_odd)

slit[indgen(nslit)].ref_ledge=lledge
slit[indgen(nslit)].ref_uedge=uuedge

for i=0,nslit-1 do begin 
 print,'slit ledge uedge =',slit[i].ref_ledge,'  ',slit[i].ref_uedge
endfor

if keyword_set(width) then begin 
 width=width 
else 
 width=1
endelse

meanx=intarr(xmax)
for i=xmax-nbs/2,nbs/2,-1 do begin 
 print,'i =',i,'  min =',i-nbs/2,' max =',i+nbs/2
 meanx[i]=i

 imystack=imy(i-nbs/2:i+nbs/2,*)
 sumstack=abs(total(imystack,1))
 

 
 for j=0,nslit-1 do begin 
   tmp_lmin=slit[j].ref_ledge-width
   tmp_lmax=slit[j].ref_ledge+width
   tmp_umin=slit[j].ref_uedge-width
   tmp_umax=slit[j].ref_uedge+width
   
   jk=max(sumstack(tmp_lmin:tmp_lmax),id)
   jk=max(sumstack(tmp_umin:tmp_umax),id2)
   slit[j].ref_ledge=idy(id+tmp_lmin)
   slit[j].ref_uedge=idy(id2+tmp_umin)
   ;print,'slit[2].ref_ledge =',slit[j].ref_ledge
   ;print,'slit[2].ref_uedge =',slit[j].ref_uedge
  ; oplot,meanx[i],slit[2].ref_ledge
 ;  oplot,meanx[i],slit[2].ref_uedge
   slit[j].ledge[i]=slit[j].ref_ledge
   slit[j].uedge[i]=slit[j].ref_uedge
 endfor

; peaksid=peaks(sumstack,200,npk)
; nbpeaks=n_elements(peaksid)
; print,'nbpeaks =',nbpeaks
; ids_eve=indgen(nslit)*2
; lledge=peaksid(ids_eve)
; ids_odd=2*indgen(nslit)+1
; uuedge=peaksid(ids_odd)
; slit[indgen(nslit)].ledge(i)=lledge
; slit[indgen(nslit)].uedge(i)=uuedge
endfor

edges_img=dblarr(2032,4064)+0.001
window,5,xsize=900,ysize=400
plot,meanx(where(slit[2].ledge gt 0)),slit[2].ledge(where(slit[2].ledge gt 0)),yrange=[940,3400]
oplot,meanx(where(slit[2].uedge gt 0)),slit[2].uedge(where(slit[2].uedge gt 0))
for i=0,nslit-1 do begin 
 oplot,meanx(where(slit[i].ledge gt 0)),slit[i].ledge(where(slit[i].ledge gt 0))
 oplot,meanx(where(slit[i].uedge gt 0)),slit[i].uedge(where(slit[i].uedge gt 0))
 edges_img(meanx(where(slit[i].uedge gt 0)),slit[i].uedge(where(slit[i].uedge gt 0)))=1.0
 edges_img(meanx(where(slit[i].ledge gt 0)),slit[i].ledge(where(slit[i].ledge gt 0)))=1.0
endfor 

; export the slit edges as fits file : 
writefits,'slit_edges.fits',edges_img

if keyword_set(norder) then begin
 norder=norder 
else 
 norder=3
endelse 

; fit 3rd-order polynomial to the slit edges along the rows : 
fit_par=dblarr(2*nslit,norder+1) 
k=0
for i=0,nslit-1 do begin
  idsok=where(slit[i].ledge gt 0)
  fit_par(k,*)=poly_fit(meanx(idsok),slit[i].ledge(idsok),norder, $ 
                         yfit=temp)
  slit[i].lfit(idsok)=temp                       
  print,'fit_par(k,*) =',fit_par(k,*)
  k+=1
  idsok=where(slit[i].uedge gt 0)
  fit_par(k,*)=poly_fit(meanx(idsok),slit[i].uedge(idsok),norder, $ 
                         yfit=temp)
   slit[i].ufit(idsok)=temp 
  k+=1
endfor

edges_fit=dblarr(xmax+1,ymax+1)+0.001
for i=0,nslit-1 do begin 
; print,'slit[i].ufit(where(slit[i].uedge gt 0)) =',slit[i].ufit(where(slit[i].uedge gt 0))
idsok=where(slit[i].uedge gt 0)
 edges_fit(meanx(idsok),floor(slit[i].ufit(idsok)))=1.0
idsok=where(slit[i].ledge gt 0)
 edges_fit(meanx(idsok),floor(slit[i].lfit(idsok)))=1.0

endfor 
writefits,'slit_edges_fit.fits',edges_fit

; writing the slitgrid : 
b_0=dblarr(2*nslit) 
b_i=dblarr(2*nslit,norder)
b_0=fit_par(*,0) 		; these are the y-intercepts  
b_i=fit_par(*,1:norder) 	; these are the other parameters of the fit. 

print,'nslit=',nslit
for i=0,2*nslit-1 do begin 
 print,'i =',i
 print,'b_i(i,*) =',b_i(i,*)
endfor
; fit b_i along b_0 (b_0 is the y-intercept) : 
q=fltarr(norder,norder+1)
for i=0,norder-1 do begin 
 q(i,*)=poly_fit(b_0,b_i(*,i),norder)
endfor  

; write the slitgrid to disk : 
x=indgen(xmax+1)
y=indgen(ymax+1)
slitarr=xy2b(x,y,q)

writefits,'slitgrid.fits',slitarr



;idymax=lclxtrem(derivy,23,/maxima)
;idymaxx=where(derivy gt 5e3)
;idd=[idymax,idymaxx]
;idd=idd(sort(idd))
;idds=shift(idd,-1)
;iddsub=idd-idds
;idok=where(iddsub eq 0) 
;edges_id=idd(idok)
;print,'edges_id before =',edges_id

; read in cross-dispersion position of the slits (roughly) : 
; the slits could also be identified using a graphical interface : 
;nslit=0
;slity=intarr(1000)
;openr,lun,'slitpos.dat',/get_lun 
; while (eof(lun) ne 1) do begin 
;   readf,lun,tmp,format='(f5.4)'
;   slity[nslit]=tmp
;   nslit=nslit+1
; endwhile 
; free_lun,lun
; ;print,'nb of slits =',nslit
; ;print,'slity =',slity(0:nslit-1);

; ; loop over each slit : 
; edges=idy(edges_id)
; slitymin=intarr(nslit) 
; slitymax=intarr(nslit)
; for i=0,nslit-1 do begin 
;    dsn=edges[where(slity[i]-edges lt 0)]
;   ; print,'dsn =',dsn 
;    subn=slity[i]-dsn
;    dsp=edges[where(slity[i]-edges gt 0)]
;  ;  print,'dsp =',dsp
;    subp=slity[i]-dsp
;    
 ;   dsna=min(abs(subn),idmin)
;    dspa=min(abs(subp),idminp) 
;    slitymin[i]=slity[i]-abs(subn(idmin))-2
;    slitymax[i]=subp(idminp)+slity[i]-2
;  ;  print,'slit ',i
 ;  ;print,'slitymin[i] =',slitymin[i]
 ;  ;print,'slitymax[i] =',slitymax[i]
   
 ;endfor

; ysize=indgen(ymax+1)+1
; slitv=intarr(ymax+1)
; for i=0,nslit-1 do begin 
;   slitv(slitymin[i]:slitymax[i])=i+1
; endfor;;

;window,2,xsize=600,ysize=400
;plot,ysize,slitv

; trace the slitedges and do a 3rd order poly fit 
; 0) create structures for nslit. 
; i) loop from x=xmax to x=0 
; ii) stack along the ydirection 
; iii) for each slit, calculate max slope in a 5 pix window centered 
;      on the value obtained previously
; iv) note the new edge position for the value of x.

; v) at the end of the calculation, you should have x vs y for all 
;    edges. Do a polynomial fit 
;nbt=10
;meanx=dblarr(2032)					; store the values of x position 
;slit=replicate({slitstrct},nslit)	; create structures for nslit ;

;for i=xmax,1000,-1 do begin
 ;print,'i =',i
 ;stack=total(imy(i-nbt:i,*),1)
 ;meanx[i]=0.5*(2.*i-nbt)
 ;print,'meanx[i] =',meanx[i]
 ;derivy=deriv(idy,stack)
 ;idymax=lclxtrem(derivy,23,/maxima)
 ;idymaxx=where(derivy gt 5e3)
 ;idd=[idymax,idymaxx]
 ;idd=idd(sort(idd))
 ;idds=shift(idd,-1)
 ;iddsub=idd-idds
 ;idok=where(iddsub eq 0) 
; edges_id=idd(idok)
; edges=idy(edges_id)
; for j=2,2 do begin
   ; print,'edges =',edges
 ;   dsn=edges[where(slity[j]-edges lt 0)]
 ;  ; print,'dsn =',dsn 
;    subn=slity[j]-dsn
;    dsp=edges[where(slity[j]-edges gt 0)]
;  ; print,'dsp =',dsp
;    subp=slity[j]-dsp
    
 ;   dsna=min(abs(subn),idmin)
 ; ;  print,'dsna =',dsna
 ;   dspa=min(abs(subp),idminp) 
 ; ;  print,'dspa =',dspa
 ;    slit[j].ledge[i]=slity[j]-dsna-2
 ; ;  print,'slit[j].ledge[i] =',slit[j].ledge[i]
 ;   slit[j].uedge[i]=dspa+slity[j]-2
 ; ;;  print,'slit[j].uedge[i] =',slit[j].uedge[i];
 ;endfor
;endfor

 ;print,'meanx =',meanx
 ;print,'slit[0].ledge =',slit[2].ledge
 ;print,'slit[0].uedge =',slit[2].uedge
; 
; window,3,xsize=900,ysize=600
; plot,meanx,slit[2].ledge,psym=2,xrange=[1000,2031],yrange=[1080,1130]
; oplot,meanx,slit[2].uedge,linestyle=1
; xsmall=meanx(where(meanx gt 1000 and meanx lt 2031)) 
; ledge_small=slit[2].ledge(where(meanx gt 1000 and meanx lt 2031))
; uedge_small=slit[2].uedge(where(meanx gt 1000 and meanx lt 2031))
 ;while 
 ;junk=poly_fit(xsmall,ledge_small,3,yfit=templ) 
 ;junk2=poly_fit(xsmall,uedge_small,3,yfit=tempu)
 ;oplot,xsmall,templ,linestyle=2
 ;oplot,xsmall,tempu,linestyle=2
 ;for i=1,nslit-1 do begin 
 ; oplot,meanx,slit[i].ledge,psym=1
 ; oplot,meanx,slit[i].uedge,psym=1
 ;endfor
 





;for i=0,n_elements(edges_id)-1 do begin 
; 
; if sumstack(edges_id(i)-1)*sumstack(edges_id(i)) lt 0 then begin 
;  edges_id2=edges_id(i)-1
; endif 
; if sumstack(edges_id(i)+1)*sumstack(edges_id(i)) lt 0 then begin 
;  edges_id2=edges_id(i)+1
; endif
;endfor
;  edges_id=[edges_id,edges_id2]
;  edges_id=edges_id(sort(edges_id))

; print,'edges_id =',edges_id

;oplot,idy(idymax),derivy(idymax),psym=4

;writefits,'junk.fits',imy



;window,1,xsize=szx,ysize=szy
;tv,imy,TRUE=interleaving[0]

;while !MOUSE.button NE 4 do begin 
; print,'Click on the first slit (at the top). Right click when done. '
; cursor,cx0,cy0,wait=4,/change,/data
; 
; ; print the cursor on the screen : 
; tmp='+'
; idcoor=strtrim(nb,2)
; xyouts,cx0-10.0,cy0-10.0,tmp,charsize=2,/device
; xyouts,cx0+15.0,cy0+15.0,idcoor,charsize=1,/device
; 
; ; trace slit : 
; wx=5   ; width in the dispersion direction. 
; xt=cx0 ; initial position where you start to trace
; 
; 
; while xt lt xmax-wx/2 do begin 
;
;   
; 
;  xt=xt+1
; endwhile 
; xt=cx0
; while xt0 gt wx/2 do begin 
; 
;   xt=xt-1
; endwhile 
; 
; 
;  
;
;  nb+=1
;endwhile 


; for all slits in the image : 

; click on traces to identify slits (start from +y to -y) 

; trace slit along the x-direction. 

; fit to the edge of the slit. 

; display the result and ask the observer if ok! 

; 






end